#ifndef METHOD_TURBO
#define METHOD_TURBO

#include "model.h"
#include "math.h"
//#include "solvers.h"
#include <iostream>
#include <omp.h>
#include <stdlib.h>
#include <mkl.h>
#include <mkl_lapacke.h>

#define PI 3.1415926535

//Cell type
enum CellType {
	Border,
	Normal
};

//Cell boundary layer solution and other information
class CellBoundarySolution
{	
public:
	int size;	//Number of points
	std::vector<int> upper_face_ind;	//array of Upper (to Bottom Face) faces indexes
	std::vector<int> upper_cell_ind;	//array of Upper (to Bottom Face) cells indexes
	std::vector<double> cell_sizes;		//array of height of cell
	int up_face;					//Up face index in array of Upper faces
	int Max_faceNumber;				//Maximum number of cell layers of structured part of the grid
	double ny_kinetic;				//average kinetic viscosity
	Face& down_face;					//Bottom face
	std::vector<double> y;	//Array of coordinats
	std::vector<double> hi;	//Array for correlations hi=-tau
	std::vector<double> U;	//Array for averaged U
	double h;				//height of turbulent layer
	double sigma;
	double U_inf;
	double U_tau;			//shear velocity
	double KARMAN;		//Karman constant
	Vector Flow_direction;	//x1 axis direction
	std::vector<double> U_av;	//array of cell average U_values
	std::vector<Vector> P_gr;	//array of P gradients
	std::map<int, var_vector> Values;	//Map from face_index to var_vector for reconstruction
	//Constructor
	CellBoundarySolution(int SolutionSize, Face& _down_face):down_face(_down_face) {
		size = SolutionSize;
		y.resize(SolutionSize);
		hi.resize(SolutionSize);
		U.resize(SolutionSize);
	};	
	//Solve model equations given external parameters
	int Solve(double c, double alpha, double gamma0, double relative_tol,double MaxIter);
	void ToTecPlot(const char *fname1, const char *fname2);	//write to .dat file
	void InitAproach();						//set hi-array for first itteration
	double norm();							//integral S(hi(y)dy) from 0 to h
	void ComputeU_1();						//function for computation U average
	void ComputeU_2();						//function for computation U when pressure X direvative equals zero
	void ComputeU_3();						//function for computation U without discontinuity on top bound for general case
	void ComputeU_4();					//law of 1/7, independening from hi
	void TrivialSolution();				//write trivial solution U = U_pr
	void ComputeU_av();					//function for computation cell averaged U_values
	double ComputeU_av_h();					//function for computation U_av value				
	void UndimU();						//calculation of dimensionless velocity; 
	void UndimU_plus();					//transition to new undimension velocity U+;
	void UndimHi();						//calculation of dimensionless correlations;
	void UndimY();						//calculation of dimensionless y;
	void UndimYplus();					//transition to new undimension cordinate y+;
	double height(double x, double Re_cr);	//calculation of height at x point for flat plane;
};

int CellBoundarySolution::Solve(double c,double alpha, double gamma0, double relative_tol, double MaxIter)
{	
	double dy = h/(size - 1);		//grid step
	for (int i=0; i<size; i++)
	{
		y[i]=i*dy;				//1D uniform grid
	};
	U_tau = sqrt(2*ny_kinetic*U_inf/h);
	double K = c*c*U_tau/(KARMAN*ny_kinetic);	
	std::vector<double> myd(size);		//turbulent viscosity
	std::vector<double> myd_dy(size);	//dirivative of turbulent viscosity
	std::vector<double> x_r(size);		//hi_r function
	for(int i=0; i<size; i++)
	{
		double gamma = gamma0 * (h*y[i] - y[i]*y[i])/(h*h);	
		double d_gamma_dy = gamma0*(h - 2*y[i])/(h*h);
		myd[i] = U_tau * h * gamma;
		myd_dy[i] = U_tau * h * d_gamma_dy;
		x_r[i] = U_tau*U_tau*gamma0*(1-(y[i]/h))/KARMAN;
	};
	
	//Banded J-Matrix
	double *ap1 = new double[size-1];	//superdiagonal elements
	double *d = new double[size];		//diagonal elements
	double *a0 = new double[size];		//constant part of diagonal elements
	double *am1 = new double[size-1];	//subdiagonal elements
	//other parameters for function LAPACKE_dgtsvx
	int* ipiv = new int[size];
	int info;
	double *dlf = new double[size-1];
	double *df = new double[size];
	double *duf = new double[size-1];
	double *du2 = new double[size-2];
	double *x = new double[size];		//solution
	double *F= new double[size];		//right hand  Jx=F or x=F/J
	double *work = new double[size * 3];	
	double rcond;
	double ferr;
	double berr;
	//calculating a0 coeficients
	a0[0] = 1;
	for(int i=1; i<(size-1); i++)
	{
		a0[i] = 2*(ny_kinetic + myd[i]) / (dy*dy);
		a0[i] = a0[i] + (ny_kinetic*alpha/(h*h) - K);
		a0[i] = a0[i] * x_r[i];
	};
	a0[size-1] = 1;
	//calculating a-1 coeficients
	for(int i=0; i<(size-2); i++)
	{
		am1[i] = myd_dy[i+1] / (2*dy);
		am1[i] = am1[i] - (ny_kinetic + myd[i+1]) / (dy*dy);
		am1[i] = am1[i] * x_r[i+1];
	};
	am1[size-2] = 0;
	//calculating a+1 coeficients
	ap1[0] = 0;
	for(int i=1; i<(size-1); i++)
	{
		ap1[i] = (-1)*myd_dy[i] / (2*dy);
		ap1[i] = ap1[i] - (ny_kinetic + myd[i]) / (dy*dy);
		ap1[i] = ap1[i] * x_r[i];
	};	
	//Newton Method
	for(int i=0; i<MaxIter; i++)
	{
		//add unsteady terms to J-matrix
		F[0] = hi[0];
		F[size-1] = hi[size-1];
		d[0] = a0[0];
		d[size-1] = a0[size-1];
		for(int j=1; j<(size-1); j++)
		{
			d[j]=a0[j] + 2*K*hi[j];
			F[j] = a0[j]*hi[j] + am1[j-1]*hi[j-1] + ap1[j]*hi[j+1] + K*hi[j]*hi[j];
		};
			
		info = LAPACKE_dgtsvx( LAPACK_COL_MAJOR,
			'N',
			'N',
			size,
			1,
			am1, d, ap1,
			dlf, df, duf,
			du2,
			ipiv,
			F, size,
			x, size,
			&rcond,
			&ferr,
			&berr
		);
		if(info!=0)
			std::cout << "U-matrix of LU decomposition is singular!LAPACKE_dgtsvx info=" << info << "\n";
		if(info<0)
		{
			std::cout << "LAPACKE_dgtsvx argument " << (-1)*info << "is illegal!\n";
			int a;
			std::cin >> a;
			break;
		};
		//relative tolerance
		double _relative_tol = 0;
		double max_ele = 0;
		//iteration cycle
		for(int j=0; j<size; j++)
		{
			hi[j] = hi[j] - x[j];
			if(_relative_tol<fabs(x[j])) _relative_tol=fabs(x[j]);
			if(max_ele<fabs(hi[j])) max_ele = hi[j];
		};
		_relative_tol = _relative_tol/max_ele;
		if(_relative_tol < relative_tol) return i;
	};
	return 0;
};
void CellBoundarySolution::ToTecPlot(const char *fname1 = "Hi_variable.dat", const char *fname2 = "U_average.dat")
{
	std::ofstream ofs(fname1);
	for(int i=0; i<size; i++)
	{
		ofs << y[i] << '\t' << hi[i] << '\n';
	};
	ofs.close();
	ofs.open(fname2);
	for(int i=1; i<size; i++)
	{
		ofs << y[i] << '\t' << U[i] << '\n';
	};
	ofs.close();
	/*ofs.open("U_correlation.dat");
	for(int i=0; i<size; i++)
	{
		ofs << y[i] << '\t' << U_correlation[i] << '\n';
	};
	ofs.close();
	ofs.open("U_Pr.dat");
	for(int i=0; i<size; i++)
	{
		ofs << y[i] << '\t' << U_pr[i] << '\n';
	};
	ofs.close();
	ofs.open("U_relative.dat");
	for(int i=0; i<size; i++)
	{
		ofs << y[i] << '\t' << U_relative[i] << '\n';
	};
	ofs.close();*/
};
void CellBoundarySolution::InitAproach()
{
	for (int i=0; i<size; i++)	
		hi[i] = sigma*h;					//initial approach 1
			//hi[i] = 0.5*sigma*h;			//initial approach 2	
			//hi[i] = sigma*h*sin(PI*((double)i)/size);
};
double CellBoundarySolution::norm()
{
	double res=0;
	double dy = h/(size - 1);
	for(int i=1; i<size; i++)
		res += 0.5*dy*(hi[i]+hi[i-1]);
	return res;
};
void CellBoundarySolution::ComputeU_1()	//TO DO CHECK
{
	//Sega part
	double tau = 0;				//integral S(tau(y))dy from 0 to y
	U[0] = 0;
	double dy = h/(size - 1);
	for(int i=1; i<size; i++)	//we accept that U[0] = 0;
	{
		tau -= 0.5*dy*(hi[i]+hi[i-1]);		//integral increment
		//correct expression for difinition U when U_inf = U_max
		//U[i] = tau/ny_kinetic - sigma*(y[i]*y[i] - 2*h*y[i])/(2*ny_kinetic);
		//correct expression for difinition U for general case
		U[i] = tau - 0.5*sigma*y[i]*y[i] + (0.5*sigma*h + ny_kinetic*U_inf/h)*y[i];
		//U[i] =  - 0.5*sigma*y[i]*y[i] + (0.5*sigma*h + ny_kinetic*U_inf/h)*y[i];
		U[i] /= ny_kinetic;
	};
	return;	
};
void CellBoundarySolution::ComputeU_2()
{
	double tau = 0;				//integral S(tau(y))dy from 0 to y
	double I = this->norm();
	I *= (-1);
	U[0] = 0;
	double dy = h/(size - 1);
	for(int i=1; i<size; i++)	//we accept that U[0] = 0;
	{
		tau -= 0.5*dy*(hi[i]+hi[i-1]);		//integral increment
		double test = ny_kinetic*U_inf;
		U[i] = tau - (I - ny_kinetic*U_inf)*y[i]/h;
		U[i] /= ny_kinetic;
	};
	return;
};
void CellBoundarySolution::ComputeU_3()
{
	////	ny_k*U(y) = S(tau(y))dy - C1*y - 0.5*sigma*y^2			////
	double tau = 0;				//integral S(tau(y))dy from 0 to y
	double I = this->norm();
	I *= (-1);
	U[0] = 0;
	double dy = h/(size - 1);
	double C1 = (I - ny_kinetic*U_inf - 0.5*sigma*h*h)/h;
	for(int i=1; i<size; i++)	//we accept that U[0] = 0;
	{
		tau -= 0.5*dy*(hi[i]+hi[i-1]);		//integral increment
		U[i] = tau - C1*y[i] - 0.5*sigma*y[i]*y[i];
		U[i] /= ny_kinetic;
	};
	return;
};
void CellBoundarySolution::ComputeU_4()
{
	double y_undim;
	for(int i=0; i<size; i++)
	{
		y_undim = y[i]/y[size-1];
		U[i] = U_inf*exp(0.143*log(y_undim));
	};
	return;
};
double CellBoundarySolution::ComputeU_av_h()
{
	double U_average = 0;
	double dy = h/size;				//grid step
	for(int i=1; i<size; i++)	//we accept that U[0] = 0;
	{
		U_average += 0.5*dy*(U[i]+U[i-1]);	
	};
	U_average /= h;
	return U_average;
};
void CellBoundarySolution::UndimU()
{
	for(int i=0; i<U.size(); i++)
	{
		U[i] /= sqrt(sigma*h);
	};
	return;
};
void CellBoundarySolution::UndimHi()
{
	for(int i=0; i<hi.size(); i++)
	{
		hi[i] /= (sigma*h);
	};
	return;
};
void CellBoundarySolution::UndimU_plus()
{
	////	tau_w/ro = -C1			////
	/*double I = this->norm();
	I *= (-1);
	double C1 = (I - ny_kinetic*U_inf - 0.5*sigma*h*h)/h;
	C1*= (-1);
	U_tau = sqrt(C1);*/
	for(int i=0; i<U.size(); i++) U[i] /= U_tau;
	return;
};
void CellBoundarySolution::UndimY()
{
	for(int i=0; i<y.size(); i++) y[i] /= h;
	return;
};
void CellBoundarySolution::UndimYplus()
{
	for(int i=0; i<y.size(); i++)
	{
		y[i] *= (U_tau/ny_kinetic);
	};
	return;
};
void CellBoundarySolution::TrivialSolution()
{
	double dy = h/(size - 1);		//grid step
	for (int i=0; i<size; i++)
	{
		y[i]=i*dy;				//1D uniform grid
		U[i] = -sigma*(y[i]*y[i] - 2*h*y[i])/(2*ny_kinetic);
		hi[i] = 0;
	};
};
void CellBoundarySolution::ComputeU_av()		//TO DO change the work of function for better precision
{
	double y_left = 0;
	int k = 0;
	double dy = y[1] - y[0];
	U_av.resize(up_face+1);
	U_av[0] = 0;
	double S_0 = 0;
	if(size<2) std::cout << "1D_grid size of turbulent equation is really low (size<2)!'\n'";
	for(int k=0, j=0; k<U_av.size(); k++)
	{
		U_av[k] += S_0;
		for(; (y[j+1]-y_left) < cell_sizes[k]; j++)
		{
			U_av[k] += 0.5*(U[j] + U[j+1])*(y[j+1] - y[j]);	//
		};
		y_left += cell_sizes[k];
		double side = U[j] + (U[j+1] - U[j])*(y_left - y[j])/(y[j+1] - y[j]);
		U_av[k] += 0.5*(side + U[j])*(y_left - y[j]);
		U_av[k] /= cell_sizes[k];
		double S_0 = 0.5*(side + U[j+1])*(y[j+1]-y_left);
		j++;
	};	
	
	U_av[k] *= (dy/cell_sizes[k]);
	double l1 = cell_sizes[k]-(y[size-1]-y_left);
	double L = l1;
	//calc distance for y[size-1] to center of up_face
	for(int i = k+1; i<=up_face; i++)
	{
		L += cell_sizes[k];
	};
	//approximation of U to top of the during cell
	double U_right = U[size-1] + (U_inf - U[size-1])*l1/L;
	U_av[k] += 0.5*(U[size-1]+U_right)*l1/cell_sizes[k];	//for linear interpolation between U_av[size-1] and U_inf
	//U_av[k] += U[size-1]*l1/cell_sizes[k];		//if U_av[size-1] is the reconstruction of U to Up_face
	double delta = (U_right - U[size-1])/l1;
	for(int i = k+1; i<=up_face; i++)
	{
		double U_left = U_right;
		U_right = U_left + delta*cell_sizes[i];
		U_av[i] = 0.5*(U_left+U_right);
	};
};
double CellBoundarySolution::height(double x, double Re_cr)	//calculation of height at x point for flat plane;
{
	if(x==0) return 0;
	double h_t, h_l, height;
	double x_cr = ny_kinetic*Re_cr/U_inf;		//point of transition to turbulent flow
	double delta_cr = 4.91*x_cr/sqrt(Re_cr);
	double delta_tur = 2*delta_cr;
	double x_tr = delta_tur*exp(0.2*log(U_inf/ny_kinetic))/0.37;
	x_tr = exp(1.25*log(x_tr));
	double x_tur = x_cr + x_tr;
	double Re = x*U_inf/ny_kinetic;
	double Re_tur = (x - x_cr)*U_inf/ny_kinetic;
	h_l = 4.9*x/sqrt(Re);
	if(x_tur<x) h_t = 0.37*(x-x_cr)/exp(0.2 * log(Re_tur));		//h_t for turbulence
	if(x_tur>=x) h_t = delta_cr + 0.5*0.37*(x-x_cr)/exp(0.2 * log(Re_tur));	//h_t for transition zone
	if(x_cr<x) height = h_t;
	else height = h_l;
	return height;
};

//Describes numerical method
class Method {
public:
	Model& model;	//Link to actual model
	Parallel& P;	//Link to wrapper for parallel part
	//Links fast access to data
	std::vector<Face>& faces;	//Global array of faces
	std::vector<Cell>& cells;	//Global array of cells
	std::vector<var_vector>& values;	//Global array of values to each cell
	int nv;	//Number of variables
	int nf; //Number of faces
	int nc; //Number of cells
	//Map indexes
	inline int IndexOf(Cell& c) { return model.grid.cell_index[c.GlobalIndex]; };
	inline int IndexOf(Face& f) { return model.grid.face_index[f.GlobalIndex]; };
	inline int IndexOf(Node& n) { return model.grid.node_index[n.GlobalIndex]; };
	inline int CellInd(int global_index) { return model.grid.cell_index[global_index]; };
	inline int FaceInd(int global_index) { return model.grid.face_index[global_index]; };	
	//Constructor
	Method(Model& _m): model(_m),
		faces(model.grid.faces),
		cells(model.grid.cells),
		values(model.values),
		nv(model.nv),
		P(model.P),
		ConfigParameters(model.ConfigParameters)
	{					
		nf = faces.size();
		nc = cells.size();
	};
	//Method parameters
	double ResidualTolerance;	//Residual tolerance
	bool SecondOrder;	//Is second order reconstruction
	int Max_Iter;	//Maximum iterations number
	double t_max;	//Calculation physical time
	double MaxExecutionTime;	//Maximum calculation time in minutes
	int OutputStep;	//Output data every N steps
	int SuppressOutput; //Is output suppressed	
	param_list& ConfigParameters;	//Configuration parameters
	int LoadConfig();	//Load config
	int Run();	//Computational cycle
	double Step(double time);	//One time step ; return time step size
	//Parallel operations
	int InitParallel(); //Init all parallel operations for method
	int ExchangeGradValues();	//Exchange gradient values
	int ExchangeLimiterValues();//Exchsnge limiter values
	//Needed variables
	std::vector<Vector> GetGradient(int Global_index); //Accessor to gradient values
	var_vector GetLimiter(int GlobalIndex);	// to limiter values
	std::map<int, std::vector<Vector> > outgrad; //Gradient values from another processes
	std::map<int, var_vector> outlimit;			//limiter values from another proesses
	std::vector<std::vector<Vector> > grad;	//Gradient in each cell
	std::vector<var_vector> limit;			//Limiters in each cell
	std::vector<var_vector> flux;	//Flux on each surface
	std::vector<var_vector> residual;	//Residual
	//Useful functions
	var_vector Limiter(Cell& c);
	std::pair<var_vector, var_vector> Method::Reconstruct(Face& f); //Reconstruction
	var_vector Reconstruct2(Cell& c, Face& f); //Complex version of previous
	var_vector InterpolationVor(Cell& c, Face& f); //Interpolate values to face for Voronoy grid
	var_vector Interpolation(Face& f);		//Interpolale values to face for any grid
	std::vector<Vector> InterGrad(Face& f);//Interpolate Gradients to face for any grid
	var_vector F(var_vector U, Vector n); //Symmetric part of flux
	std::pair<var_vector, double> Flux(Face& f); //Solve Reimann problem
	var_vector VolumeSource(Cell& c);	//Volume sources
	var_vector ViscosityFlux(Face& f); //Calculate diffusion part of flux
	double Harten(double z, double eps); //fabs() and Harten's entropy correction procedure
	std::vector<Vector> Grad(Cell& c); //Calculate solution gradient for given cell
	std::vector<Cell> Closest(Cell& c, int n);	//Find closest neighbours of given cell
	//Boundary layer treatment functions and variables
	bool isBoundaryLayerActive;		//Is boundary layer modeled
	double StartTurbulenceTime;		//Time when turbulence may start
	double BoundaryAlpha;			//Parameters
	double BoundaryC;				//
	double BoundaryGamma;			//
	int MaxBoundaryCells;			//
	double BoundaryMaxIter;			//maximum number of iterations for correlations computation 
	double BoundaryRelativeTol;		//relative tolerance of solution
	std::vector<int> BoundaryMarkers;	//Markers for walls
	//std::vector<CellType> cell_types;	//Cell type
	//std::vector<double> cell_distance;	//Distance to closest wall
	std::map<int, CellBoundarySolution*> solutions; //Hash of solutions for boundary layer Face_index -> solution
	void InitBoundaryLayer();			//Init boundary layer model parameters
	void SolveBoundaryLayer(Cell& c, Vector Uinf);	//Solve boundary layer model for given cell
	int GetNextFace(CellBoundarySolution& our_cell);		//define the upper face if it is't known for orthogonal grid
	inline int GetPreviousFace(CellBoundarySolution& our_cell);//increment of up_face index
	double DefineUinf(CellBoundarySolution& our_cell);		//define U_inf
	void CalcP_grad(CellBoundarySolution& our_cell, int Cell_number);		//calc & write P_gradients from bottom to Face_number cell
	Vector ChooseFlow_direction(CellBoundarySolution& our_cell);//define Flow Direction for turbulent flow
																//find average sigma and average kinetic viscocity ny_kinetic
	bool Check(CellBoundarySolution& our_cell);			//check that all parameters are well prepared
	bool CheckTurbulent(CellBoundarySolution& our_cell);	//check is flow turbulent or not
	
	//Volume sources treatment	
	bool isVolumeSourcesActive; //Is volume sources active	
	std::map<std::string, bool> VolumeSourcesTypes; //List of source types currently active
	void InitVolumeSources(); //Read config and load parameters	
	var_vector CalcVolumeSources(int GlobalIndex); //Calculace volume sources term for given cell
	//List of parameters for volume sources	
	Vector g; //Uniform gravity	vector	
};

///Boundary layer treatment functions
//TO DO create opposite function
int Method::GetNextFace(CellBoundarySolution& our_cell)
{
	int upperSize = our_cell.upper_face_ind.size(); //size of array of already calculated faces indexes
	//GetNextFace calls first time - initialisation of up_face
	if(upperSize==0) our_cell.up_face = -1;
	//the case of pass out of structured part of the grid
	if(our_cell.up_face==(our_cell.Max_faceNumber-1))
	{
		std::cout << "Maximum number of layers has been reached at face " << our_cell.down_face.GlobalIndex << '\n';
		int read;
		std::cin >> read;
		exit(0);
	};
	//the case when new face is already calculated	
	if ((our_cell.up_face!=-1)&&(our_cell.up_face < (upperSize-1))) {
		our_cell.up_face++;
		return our_cell.upper_face_ind[our_cell.up_face];
	};
	// find new up face
	//TO DO do generalize for MPI version
	//TO DO implement Model::GetFace()
	Face f;				//Up face 
	int cell_up_upper;	//GlobalIndex of cell between up and upper faces
	// up_face = the bottom face
	if(our_cell.up_face==-1)
	{
		f = our_cell.down_face;
		cell_up_upper = f.FaceCell_1;
	} else {
		f = faces[FaceInd(our_cell.upper_face_ind[our_cell.up_face])];	//old up_face
		//GlobalIndex of cell placed under up_face
		int cell_up_down = our_cell.upper_cell_ind[our_cell.up_face];
		if(f.FaceCell_1==cell_up_down) cell_up_upper = f.FaceCell_2;
		else{
			cell_up_upper = f.FaceCell_1;
		};
	};
	//Cell between up and upper faces
	Cell c = cells[CellInd(cell_up_upper)];
	//dimension of the problem
	int Dim;
	if(f.FaceNodes.size()==2) Dim = 2;	//for 2D case
	if(f.FaceNodes.size()>2) Dim=3;		//for 3D case
	if(f.FaceNodes.size()<2) std::cout << "Warning: face has less the 2 nodes!" << '\n';
	std::vector<int> upper_nodes;	//3 or 2 nodes the upper_face includes
	for(int i=0, j; upper_nodes.size()<Dim;i++)
	{
		for( j=0; j<f.FaceNodes.size();j++) if(c.Nodes[i]==f.FaceNodes[j]) break;
		if(j==f.FaceNodes.size()) upper_nodes.push_back(c.Nodes[i]);
	};
	Face upper_face;
	//search face that containes all upper_nodes
	for(int i=0, num=0; num<Dim; i++)
	{
		num = 0;
		upper_face = faces[FaceInd(c.Faces[i])];
		for(int j=0, k=0; k<upper_face.FaceNodes.size(); j++)
		{
			for(k=0; k<upper_face.FaceNodes.size(); k++) 
				if(upper_face.FaceNodes[k]==upper_nodes[j]) 
				{
					num++;
					break;
				};
			if(num==Dim) break;
		};
	};
	//It's may be good to generalize the case below for when y_direction can be significantly changed from one to other up_face
	Vector y_direction = (upper_face.FaceCenter - our_cell.down_face.FaceCenter);	//y_direction
	y_direction /= y_direction.mod();	//unit vector of direction
	double segment = ((upper_face.FaceCenter - f.FaceCenter)*y_direction);		//projection of new segment to y_direction
	our_cell.upper_face_ind.push_back(upper_face.GlobalIndex);
	our_cell.upper_cell_ind.push_back(cell_up_upper);
	our_cell.cell_sizes.push_back(segment);
	our_cell.up_face++;
	return upper_face.GlobalIndex;
};
inline int Method::GetPreviousFace(CellBoundarySolution& our_cell)
{
	if(our_cell.up_face==0)
	{
		std::cout << "Can't GetPreviousFace\n" ;
	};
	return our_cell.up_face--;
};
double Method::DefineUinf(CellBoundarySolution& our_cell)
{
	var_vector Values(model.nv);		//values from cell placed upper top face - Down_Cell
	Face up_face = faces[FaceInd(our_cell.upper_face_ind[our_cell.up_face])];
	int down_cell = our_cell.upper_cell_ind[our_cell.up_face];	//index of Down_Cell
	if(up_face.FaceCell_1==down_cell) Values = Reconstruct(up_face).second;
	else Values = Reconstruct(up_face).first;
	Vector Velocity_inf(Values[1],Values[2],Values[3]);	//ro*V
	Velocity_inf /= Values[0];							//infinum velocity
	Velocity_inf = (Velocity_inf*our_cell.Flow_direction)*our_cell.Flow_direction;	//projection of infinum Velocity to flow direction
	return Velocity_inf.mod();
};
//calculate grads of pressure for cells array from bottom to Max_faceNumber cell
void Method::CalcP_grad(CellBoundarySolution& our_cell, int Cell_number)
{
	our_cell.P_gr.resize(Cell_number);	//array of gradients(P)
	std::vector<Vector> gradCell;								//gradients in cell[i]
	for(int i=0; i<our_cell.P_gr.size(); i++)		
	{
		gradCell = grad[CellInd(our_cell.upper_cell_ind[i])];
		our_cell.P_gr[i] = gradCell[4];	// grad(ro*E)
		Vector Velocity(values[CellInd(our_cell.upper_cell_ind[i])][1],values[CellInd(our_cell.upper_cell_ind[i])][2],values[CellInd(our_cell.upper_cell_ind[i])][3]);	//ro*velocity in cell[i]
		Velocity /= values[CellInd(our_cell.upper_cell_ind[i])][0];		//velocity in cell[i]		ro*grad(V^2) = ro*u*grad_u + ro*v*grad_v + ro*w*grad_w
		our_cell.P_gr[i] -= Velocity.x*(gradCell[1]-Velocity.x*gradCell[0]);			//grad(ro*E) - ro*u*grad_u				
		our_cell.P_gr[i] -= Velocity.y*(gradCell[2]-Velocity.y*gradCell[0]);			//grad(ro*E) - ro*u*grad_u - ro*v*grad_v
		our_cell.P_gr[i] -= Velocity.z*(gradCell[3]-Velocity.z*gradCell[0]);			//grad(ro*E) - ro*u*grad_u - ro*v*grad_v - ro*w*grad_w
		our_cell.P_gr[i] -= 0.5*(Velocity*Velocity)*grad[model.grid.cell_index[our_cell.upper_cell_ind[i]]][0];	 // grad(ro*E) - 0.5*V^2*ro - ro*u*grad_u - ro*v*grad_v - ro*w*grad_w
		our_cell.P_gr[i] *= model.gamma-1;					//gradient of P for cell[i]
	};
};
////TO DO CHECK, THE FUNCTION RETORN ZERO VECTOR
Vector Method::ChooseFlow_direction(CellBoundarySolution& our_cell)
{
	Vector res(0,0,0);
	Vector down_point = our_cell.down_face.FaceCenter;		//ends of segment form one to another face
	double ro_ave = 0;				//average density
	Vector up_point = faces[FaceInd(our_cell.upper_face_ind[our_cell.up_face])].FaceCenter;
	Vector y_direction = (up_point - our_cell.down_face.FaceCenter);		//direction of y axy
	double Length = y_direction.mod();				//y_direction as unity vector
	y_direction /= Length;
	for(int i=0; i<=our_cell.up_face; i++)
	{
		res += our_cell.P_gr[i]*our_cell.cell_sizes[i];			//P_average gradient = P grad[i] * l[i]
		ro_ave += values[CellInd(our_cell.upper_cell_ind[i])][0]*our_cell.cell_sizes[i];	//Ro[i]*l[i]
	};
	res -= (res*y_direction)*y_direction;	//flow direction is perpendicalar to y_direction
	our_cell.sigma = res.mod()/ro_ave;	//sigma = (grad_P)_average/ro_average
	ro_ave /= Length;		//really average density
	our_cell.ny_kinetic = model.my/ro_ave;	//average kinetic viscosity
	if(res.mod()!=0) res /= res.mod();		//flow_direction is a unit vector
	return -res;
};
//check the condition that turbulent layer height is down than our choosed height 
bool Method::Check(CellBoundarySolution& our_cell)
{
	if(our_cell.sigma==0)	//if no pressure gradient
	{	
		our_cell.h=0;
		return true;
	};
	our_cell.h = sqrt(2*our_cell.ny_kinetic*our_cell.U_inf/our_cell.sigma);
	double length = ((faces[FaceInd(our_cell.upper_face_ind[our_cell.up_face])]).FaceCenter - our_cell.down_face.FaceCenter).mod();
	if(our_cell.h<=length) return true;
	return false;
};

//check that flow is laminar
bool Method::CheckTurbulent(CellBoundarySolution& our_cell)
{
	if(our_cell.h==0) {
		return false;
	}
	else
	{
		double Re_our = our_cell.h*sqrt(our_cell.h*our_cell.sigma)/our_cell.ny_kinetic;
		if(Re_our<sqrt(BoundaryAlpha/(BoundaryC*BoundaryC))) return false;	//we take Re critical as alfa/c^2
		else return true;
	};
};

//Get gradient value
std::vector<Vector> Method::GetGradient(int Global_index) {
	if (P.IsLocal(Global_index)) return grad[model.grid.cell_index[Global_index]];
	if (outgrad.find(Global_index) != outgrad.end()) {		
		return outgrad[Global_index];
	};
	throw 1;	
	std::vector<Vector> res;
	return res;
};

//Get limiter value
var_vector Method::GetLimiter(int Global_index){
	if (P.IsLocal(Global_index)) return limit[model.grid.cell_index[Global_index]];
	if (outlimit.find(Global_index) != outlimit.end()) {		
		return outlimit[Global_index];
	};
	throw 1;	
	var_vector res;
	return res;
};

//Gradient calculation for a given cell
std::vector<Vector> Method::Grad(Cell& c)
{		
	std::vector<Vector> grad(nv);
	const int m_max = 10;	
	const int nrhs_max = 6;
	//For all values write linear system like A*grad = B;	
	MKL_INT m = c.Faces.size();
	const MKL_INT n = 3;
	const MKL_INT nrhs = nv;
	double A[m_max*n];	
	int lda = m;
	//write matrix constant A depending only on our grid
	for(int i=0; i<m; i++)
	{
		Face cur_face = faces[FaceInd(c.Faces[i])];
		//for boundary face
		if(cur_face.isExternal)
		{
			A[i] = (cur_face.FaceCenter - c.CellCenter).x;
			A[i+lda] = (cur_face.FaceCenter - c.CellCenter).y;
			A[i+2*lda] = (cur_face.FaceCenter - c.CellCenter).z;
			continue;
		};
		Cell c_neigh;
		if (cur_face.FaceCell_1 == c.GlobalIndex)
		{
			c_neigh = model.GetCell(cur_face.FaceCell_2);
		}	else {
			c_neigh = model.GetCell(cur_face.FaceCell_1);
		};
		A[i] = (c_neigh.CellCenter - c.CellCenter).x;
		A[i+lda] = (c_neigh.CellCenter - c.CellCenter).y;
		A[i+2*lda] = (c_neigh.CellCenter - c.CellCenter).z;			
	};	
	double B[m_max*nrhs_max];
	int ldb = m;
	//definition of right matrix B 
	for(int i=0; i<nrhs; i++)
	{
		for(int j = 0; j<m; j++)
		{
			Face cur_face = faces[FaceInd(c.Faces[j])];
			//for boundary face
			if(cur_face.isExternal == 1)
			{
				B[j+ldb*i] = model.GetBoundaryVal(cur_face, model.GetValue(c.GlobalIndex))[i] - model.GetValue(c.GlobalIndex)[i];
				continue;
			};
			Cell c_neigh;
			if (cur_face.FaceCell_1 == c.GlobalIndex)
			{
				c_neigh = model.GetCell(cur_face.FaceCell_2);
			}else c_neigh = model.GetCell(cur_face.FaceCell_1);
			B[j+ldb*i] = model.GetValue(c_neigh.GlobalIndex)[i] - model.GetValue(c.GlobalIndex)[i];
		};				
	}	
	const int lwork = 1000;
	double work[lwork];
	MKL_INT info;
	//leasted squares solver
	//print_matrix( "A", m, n, A, lda );
	//print_matrix( "B", m, nrhs, B, ldb );	
	int rank;
	double rcond = -1.0;
	double s[m_max];
	info = LAPACKE_dgelsd( LAPACK_COL_MAJOR, m, n, nrhs, A, lda, B, ldb, s, rcond, &rank );
	//print_matrix( "Least squares solution", n, nrhs, B, ldb );
	//printf("info = %ld\n", info);
	//if (c.GlobalIndex == 200) exit(0); 
	
	if (info != 0) {		
		throw Exception("Least squares failed");
	};

	//write solution
	for(int i=0; i<nv; i++)
	{
		grad[i].x = B[ldb*i];
		grad[i].y = B[ldb*i + 1];
		grad[i].z = B[ldb*i + 2];
	};	
	return grad;	
}

//Init parallel operations specific to method
int Method::InitParallel() {
	P.AddRequestGroup("grad");
	P.AddRequestGroup("limit");
	std::vector<std::pair<int, int> > cells_to_send(0);	//what cell to send and where to send	
	std::vector<std::pair<int, int> > cells_to_recv(0); //what cell to recv and from whom to reciev		
	for (int i = 0; i<model.grid.faces.size(); i++) {
		Face Fi = model.grid.faces[i];
		if (Fi.isExternal) continue;	//Ignore external faces
		bool isFace1 = P.IsLocal(Fi.FaceCell_1);
		bool isFace2 = P.IsLocal(Fi.FaceCell_2);		
		if ((isFace1) && (!isFace2)) {
			std::pair<int, int> cell_dest;
			cell_dest.first = Fi.FaceCell_1;
			cell_dest.second = P.GetRank(Fi.FaceCell_2);
			cells_to_send.push_back(cell_dest);

			std::pair<int, int> cell_recv;
			cell_recv.first = Fi.FaceCell_2;
			cell_recv.second = P.GetRank(Fi.FaceCell_2);			
			cells_to_recv.push_back(cell_recv);
		};
		if ((!isFace1) && (isFace2)) {
			std::pair<int, int> cell_dest;
			cell_dest.first = Fi.FaceCell_2;
			cell_dest.second = P.GetRank(Fi.FaceCell_1);
			cells_to_send.push_back(cell_dest);			
			
			std::pair<int, int> cell_recv;
			cell_recv.first = Fi.FaceCell_1;
			cell_recv.second = P.GetRank(Fi.FaceCell_1);			
			cells_to_recv.push_back(cell_recv);
		};
	};	
	//Initialize send requests
	for (int i  = 0; i<cells_to_send.size(); i++) {		
		P.request_groups["grad"]->AddSendRequest(cells_to_send[i].first, nv*3, cells_to_send[i].second, cells_to_send[i].first);	
		P.request_groups["limit"]->AddSendRequest(cells_to_send[i].first, nv, cells_to_send[i].second, cells_to_send[i].first);	////???
	};		
	//Initialize recv requests
	for (int i  = 0; i<cells_to_recv.size(); i++) {	
		P.request_groups["grad"]->AddRecvRequest(cells_to_recv[i].first, nv*3, cells_to_recv[i].second, cells_to_recv[i].first);	
		P.request_groups["limit"]->AddRecvRequest(cells_to_recv[i].first, nv, cells_to_recv[i].second, cells_to_recv[i].first);////???
	};	
	outgrad.clear();	
	return 0;
};


//Exchange gradient values
int Method::ExchangeGradValues() {

	//First we send all we need	
	for (int i = 0; i<P.request_groups["grad"]->requests_to_send.size(); i++) {		
		std::vector<Vector> val = GetGradient(P.request_groups["grad"]->requests_to_send[i]->cell_index);
		for (int j = 0; j<val.size(); j++) {
			P.request_groups["grad"]->requests_to_send[i]->buf[j*3] = val[j].x;
			P.request_groups["grad"]->requests_to_send[i]->buf[j*3+1] = val[j].y;
			P.request_groups["grad"]->requests_to_send[i]->buf[j*3+2] = val[j].z;
		};
	};	
	P.request_groups["grad"]->SendAll();
	//Then we recv all needed data
	P.request_groups["grad"]->RecvAll();

	outgrad.clear();
	for (int i = 0; i<P.request_groups["grad"]->requests_to_recv.size(); i++) {		
		std::vector<Vector>  val(nv);
		for (int j = 0; j<val.size(); j++) {
			val[j].x = P.request_groups["grad"]->requests_to_recv[i]->buf[j*3];
			val[j].y = P.request_groups["grad"]->requests_to_recv[i]->buf[j*3+1];
			val[j].z = P.request_groups["grad"]->requests_to_recv[i]->buf[j*3+2];
		};		
		outgrad[P.request_groups["grad"]->requests_to_recv[i]->cell_index] = val;
	};
	//Wait for others to complete
	P.Barrier();
	return 0;
};
int Method::ExchangeLimiterValues() {

	//First we send all we need	
	for (int i = 0; i<P.request_groups["limit"]->requests_to_send.size(); i++) {		
		var_vector val = GetLimiter(P.request_groups["limit"]->requests_to_send[i]->cell_index);
		for (int j = 0; j<val.size(); j++) {
			P.request_groups["limit"]->requests_to_send[i]->buf[j] = val[j];
		};
	};	
	P.request_groups["limit"]->SendAll();
	//Then we recv all needed data
	P.request_groups["limit"]->RecvAll();

	outlimit.clear();
	for (int i = 0; i<P.request_groups["limit"]->requests_to_recv.size(); i++) {		
		var_vector val(nv);
		for (int j = 0; j<val.size(); j++) {
			val[j] = P.request_groups["limit"]->requests_to_recv[i]->buf[j];
		};		
		outlimit[P.request_groups["limit"]->requests_to_recv[i]->cell_index] = val;
	};
	//Wait for others to complete
	P.Barrier();
	return 0;
};

//Limiter function
var_vector Method::Limiter(Cell& c)
{
	var_vector U_c(nv);									//Values of c cell
	U_c = model.GetValue(c.GlobalIndex);
	var_vector U_max,U_min;								//Max and Min Values from c cell and its neighbours
	U_max = U_c;
	U_min = U_c;
	for (int i = 0; i<c.Faces.size(); i++)				//find U_min and U_max
	{
		var_vector U_neigh(nv);							//Values in i-neighbour cell
		Cell neigh;
		Face cur_face = faces[FaceInd(c.Faces[i])];
		if(cur_face.isExternal!=0) continue;
		if(cur_face.FaceCell_1 == c.GlobalIndex)		//find i-neighbour cell
		{
			neigh = cells[CellInd(cur_face.FaceCell_2)];
		}else neigh = cells[CellInd(cur_face.FaceCell_1)];
		U_neigh = model.GetValue(neigh.GlobalIndex);	//Values of i-neighbour cell
		for (int j = 0; j < nv; j++)
		{
				if(U_max[j] < U_neigh[j]) U_max[j] = U_neigh[j];	//find U_max,U_min
				if(U_min[j] > U_neigh[j]) U_min[j] = U_neigh[j];
		};
	};
	var_vector res(nv,1);						
	std::vector<Vector> grad(nv);						//Gradient in c cell
	grad = GetGradient(c.GlobalIndex);
	for(int i=0; i<c.Faces.size(); i++)				//calc limiter - res
	{
		var_vector res_new(nv);
		Face cur_face = faces[FaceInd(c.Faces[i])];
		var_vector delta2(nv);							//special function					
		for(int j=0; j<nv; j++)
		{
			delta2[j] = grad[j]*(cur_face.FaceCenter - c.CellCenter);
			if(delta2[j] > 0) res_new[j] = (U_max[j] - U_c[j])/delta2[j];
			if(delta2[j] < 0) res_new[j] = (U_min[j] - U_c[j])/delta2[j];
			if(delta2[j] ==0) res_new[j] = 1;
			if(res_new[j] > 1)	res_new[j] = 1;
			if(res[j] > res_new[j]) res[j] = res_new[j];	//find min res_new
		};
	};
	return res;
};

//Reconstruct right and left values on given face
std::pair<var_vector, var_vector> Method::Reconstruct(Face& f)
{	
	var_vector UL;
	var_vector UR;
	if (!f.isExternal) {
		UL = model.GetValue(f.FaceCell_1);
		UR = model.GetValue(f.FaceCell_2);	
	} else {
		UL = model.GetValue(f.FaceCell_1);
		UR = model.GetBoundaryValue(f, UL);
	};		
	//First order reconstruction
	if (!SecondOrder) {
		return std::pair<var_vector, var_vector>(UL, UR);		
	};
	//Second order reconstruction
	if (!f.isExternal) {
		std::vector<Vector> GradL = GetGradient(f.FaceCell_1);
		std::vector<Vector> GradR = GetGradient(f.FaceCell_2);
		Vector drL = f.FaceCenter - model.GetCell(f.FaceCell_1).CellCenter;
		Vector drR = f.FaceCenter - model.GetCell(f.FaceCell_2).CellCenter;
		var_vector dUL(nv);
		var_vector dUR(nv);
		for (int i = 0; i<dUL.size(); i++) {
			dUL[i] = drL * GradL[i];
			dUL[i] = dUL[i]*GetLimiter(f.FaceCell_1)[i];
			dUR[i] = drR * GradR[i];
			dUR[i] = dUR[i]*GetLimiter(f.FaceCell_2)[i];
		};
		UL += dUL;
		UR += dUR;
	}else
	{
		std::vector<Vector> GradL = GetGradient(f.FaceCell_1);
		Vector drL = f.FaceCenter - model.GetCell(f.FaceCell_1).CellCenter;
		var_vector dUL(nv);
		for (int i = 0; i<dUL.size(); i++) {
			dUL[i] = drL * GradL[i];
			dUL[i] = dUL[i]*GetLimiter(f.FaceCell_1)[i];
		};
		UL += dUL;
	};
	return std::pair<var_vector, var_vector>(UL, UR);
};

//Helpful functions
//Flux through cell surface
var_vector Method::F(var_vector U, Vector n)
{
	double ro = U[0];
	double vx = U[1]/ro;
	double vy = U[2]/ro;
	double vz = U[3]/ro;
	double roe = U[4];	//ro*e
	double p = (model.gamma-1.0)*(roe-ro*(vx*vx+vy*vy+vz*vz)/2.0);
	var_vector res(nv,0);
	res[0] = n.x * (ro*vx) + n.y*(ro*vy) + n.z*(ro*vz);
	double vn = vx*n.x + vy*n.y + vz*n.z;
	res[1] = n.x * (ro*vx*vx+p) + n.y*(ro*vx*vy) + n.z*(ro*vx*vz);
	res[2] = n.x * (ro*vy*vx) + n.y*(ro*vy*vy+p) + n.z*(ro*vy*vz);
	res[3] = n.x * (ro*vz*vx) + n.y*(ro*vz*vy) + n.z*(ro*vz*vz+p);
	res[4] = (n.x * vx + n.y * vy + n.z * vz)*(roe+p);
	//for (int i = 0; i<nv; i++) printf("1:%lg\n", res[i]);
	return res;
};

//fabs() and Harten's entropy correction procedure
double Method::Harten(double z, double eps) 
{
	z = fabs(z);
	if (z<eps) z = ((z*z)/eps + eps)*0.5;
	return z;
};

const double bad_volume = 0;

//Compute flux through given face
//Returns also maximal eigenvalue of Jacobian matrix
std::pair<var_vector, double> Method::Flux(Face& f)
{
	var_vector res(nv, 0);	
	//Reconstruct left and right solution parts at edge	
	//Cell c_l = cells[CellInd(f.FaceCell_1)];	//Left cell
	//var_vector UL = Reconstruct(c_l, f);
	//TO DO no reconstruction atm
	std::pair<var_vector, var_vector> rec = Reconstruct(f);
	var_vector UL = rec.first;
	var_vector UR = rec.second;
	//Calculate symmetric flux part
	//for (int i = 0; i<nv; i++) printf("%lg\n", UL[i]);
	//res = F(UL, f.FaceNormal);
	//for (int i = 0; i<nv; i++) printf("%lg\n", res[i]);
	res += 1.0*(F(UL, (1.0)*f.FaceNormal) + F(UR, (1.0)*f.FaceNormal)); //is 1.0 or 0.5 TODO		
	//for (int i = 0; i<nv; i++) printf("%lg\n", res[i]);
	//Calculates stabilization term which is a part of numerical
	//flux vector i.e. |A|(Q{R}-Q{L})
	// Roe type averaging procedure first
	double ro_l = UL[0];
	double ro_r = UR[0];
	Vector velocity_l, velocity_r;
	velocity_l.x = UL[1]/ro_l;
	velocity_l.y = UL[2]/ro_l;
	velocity_l.z = UL[3]/ro_l;
	velocity_r.x = UR[1]/ro_r;
	velocity_r.y = UR[2]/ro_r;
	velocity_r.z = UR[3]/ro_r;
	double e_l = UL[4]/ro_l;
	double e_r = UR[4]/ro_r;
	double k;
	k = 0.5*(velocity_l*velocity_l);		//kinetik energy
	double h_l = (e_l-k)*model.gamma + k;	//enthalpy
	k = 0.5*(velocity_r*velocity_r);		//kinetik energy
	double h_r = (e_r-k)*model.gamma + k;	//enthalpy
	double ro = sqrt(ro_l*ro_r);          // (Roe averaged) density
	//Debug
	for (int k = 0; k<nv; k++) if (res[k] != res[k]) {
		printf("!!!\n");
		exit(0);
	};
	double ql  = sqrt(ro_l)/(sqrt(ro_l)+sqrt(ro_r));
	double qr  = sqrt(ro_r)/(sqrt(ro_l)+sqrt(ro_r));
	Vector velocity = ql*velocity_l + qr*velocity_r;	// (Roe averaged) velocity	
	double h  = ql*h_l + qr*h_r;  // (Roe averaged) total enthalpy
	//Proceed to solution
	double phi2 = 0.5*(model.gamma - 1)*(velocity*velocity);
	double dn = f.FaceNormal.mod();
	double c = sqrt((model.gamma - 1)*h - phi2);	//acoustic velocity
	//Debug	
	double uw = velocity * f.FaceNormal;
	double eig_max = fabs(uw)+c*dn;	
    double AA1 = Harten(uw, model.eps*eig_max);       // AA1, AA3, AA1 -
    double AA3 = Harten(uw+c*dn, model.eps*eig_max);  // eigenvalues of a flux vector
	double AA4 = Harten(uw-c*dn, model.eps*eig_max);  // Jacobian matrix
    double Eig1= AA1;
    double Eiga=(AA3 - AA4)*0.5/(dn*c);
    double Eigb=(AA3 + AA4)*0.5 - AA1;
	//parametrs vectors Qa and Qb (i guess)
	var_vector Qa(nv, 0);
	var_vector Qb(nv, 0);
	Qa[0]=0;
	Qa[1]=f.FaceNormal.x;
	Qa[2]=f.FaceNormal.y;
	Qa[3]=f.FaceNormal.z;
    Qa[4]=uw;    
    Qb[0]=1;
	Qb[1]=velocity.x;
	Qb[2]=velocity.y;
	Qb[3]=velocity.z;
    Qb[4]=h;
	//Calculate solution
	//Some quotients
	double R1 =phi2*ro_r-(model.gamma-1)*ro_r*(velocity*velocity_r-e_r);	//PR =R1
	double D1 =ro_r*f.FaceNormal*(velocity_r - velocity);				//DR =D1

	double R2 =phi2*ro_l-(model.gamma-1)*ro_l*(velocity*velocity_l-e_l);	//PR =R1
	double D2 =ro_l*f.FaceNormal*(velocity_l - velocity);		

    double C2 = 1.0/(c*c);
    double DN2= 1.0/(dn*dn);
	for(int i=0; i<nv; i++){
			  res[i]+=(Eig1*UL[i] + Eiga*(R2*Qa[i]     +D2*Qb[i])
								+ Eigb*(R2*Qb[i]*C2  +D2*Qa[i]*DN2)
					-Eig1*UR[i] - Eiga*(R1*Qa[i]     +D1*Qb[i])
								- Eigb*(R1*Qb[i]*C2  +D1*Qa[i]*DN2));
			  res[i]*= 0.5;
	};	

	
	return std::pair<var_vector, double>(res, eig_max);
};

var_vector Method::InterpolationVor(Cell& c, Face& f)// MAY BE DELETED //values interpolation from cells c to inner face f
{
	var_vector inter(nv,0);
	double l_fc = (f.FaceCenter - c.CellCenter).mod();
	double l_cc = (model.GetCell(f.FaceCell_1).CellCenter - model.GetCell(f.FaceCell_2).CellCenter).mod();
	inter = model.GetValue(f.FaceCell_2) - model.GetValue(f.FaceCell_1);
	if ((f.FaceCell_1)!=c.GlobalIndex) inter=(-1)*inter;
	inter = (l_fc/l_cc)*inter;
	inter += model.GetValue(c.GlobalIndex);
	return inter;
}
var_vector Method::Interpolation(Face& f)		// values interpolation to inner face f
{
	var_vector inter(nv,0);
	Vector c_1 = model.GetCell(f.FaceCell_1).CellCenter;
	if(f.FaceCell_2==(-1)) printf("Interpolation for Values - Er");
	Vector c_2 = model.GetCell(f.FaceCell_2).CellCenter;
	double l_cc = (c_1 - c_2).mod();
	double l_fc = l_cc*(f.FaceNormal*(f.FaceCenter - c_1))/(f.FaceNormal*(c_2 - c_1));
	inter = model.GetValue(f.FaceCell_2) - model.GetValue(f.FaceCell_1);
	inter = (l_fc/l_cc)*inter;
	inter += model.GetValue(f.FaceCell_1);
	return inter;
}
std::vector<Vector> Method::InterGrad(Face& f)
	// gradients interpolation to inner face f
{
	Vector zero(0,0,0);
	std::vector<Vector> inter(nv, zero);
	Vector c_1 = model.GetCell(f.FaceCell_1).CellCenter;
	if(f.FaceCell_2==(-1)) printf("Interpolation for Grads - Er");
	Vector c_2 = model.GetCell(f.FaceCell_2).CellCenter;
	double l_cc = (c_1 - c_2).mod();
	double l_fc = l_cc*(f.FaceNormal*(f.FaceCenter - c_1))/(f.FaceNormal*(c_2 - c_1));
	for(int i=0; i<nv; i++) 
	{
		inter[i] = GetGradient(f.FaceCell_2)[i] - GetGradient(f.FaceCell_1)[i];
		inter[i] = (l_fc/l_cc)*inter[i];
		inter[i] = inter[i] + GetGradient(f.FaceCell_1)[i];
	};
	return inter;
};

var_vector Method::ViscosityFlux(Face& f)
{	
	var_vector res(nv, 0);		//result for ViscosityFlux
	std::vector<Vector> gr(nv);	//gradient on face f
	var_vector a(nv);				//for inter(extra)polation U
	
	//calculation of gradient on face f - gr
	if (!f.isExternal) {		//for inner faces
		a = Interpolation(f);
		gr =InterGrad(f);		
	} else {					//for external faces
		return res;				//TO DO
		Cell c_1;														//inner cell
		std::vector<Vector> grad1;										//inner gradient
		c_1 = model.GetCell(f.FaceCell_1);
		a = model.GetBoundaryVal(f , model.GetValue(c_1.GlobalIndex));	//values on external face f
		grad1 = GetGradient(f.FaceCell_1);	
		gr = model.GetBoundaryGradient(f, c_1, model.GetValue(c_1.GlobalIndex), grad1);
	};

	double ny = model.my;		//common viscosity = dynamic viscosity
	double Pr_com = model.Pr;	//common Prandtl number
	//double Pr_com = (ny/model.my)*model.Pr //for NS+
	Vector v_i(1,0,0);
	Vector v_j(0,1,0);
	Vector v_k(0,0,1);
	/*double C;				//divergence of velocity
	C = gr[1].x+gr[2].y+gr[3].z;
	//calc viscosity fluxes for velocity
	Vector tau_x(gr[1].x,gr[2].x,gr[3].x);
	tau_x = tau_x + gr[1];
	Vector tau_y(gr[1].y,gr[2].y,gr[3].y);
	tau_y = tau_y + gr[2];
	Vector tau_z(gr[1].z,gr[2].z,gr[3].z);
	tau_z = tau_z + gr[3];
	res[1] = ny*(((2/3)*C*v_i) - tau_x)*f.FaceNormal;
	res[2] = ny*(((2/3)*C*v_j) - tau_y)*f.FaceNormal;
	res[3] = ny*(((2/3)*C*v_k) - tau_z)*f.FaceNormal;
	//calc viscocity flux for energy
	Vector U(a[1],a[2],a[3]);
	Vector GradU2 = 2*(a[1]*gr[1] + a[2]*gr[2] + a[3]*gr[3]);
	Vector p_e = 0.5*(model.gamma/model.Pr - 1)*GradU2 -(1/3)*C*U - (model.gamma/(Pr_com))*gr[4];
	res[4] = ny*p_e*f.FaceNormal;*/
	Vector U(a[1]/a[0],a[2]/a[0],a[3]/a[0]);				// velocity on face
	double ro_C = (gr[1].x+gr[2].y+gr[3].z)-(gr[0]*U);		// ro*div(U)
	//calc viscosity fluxes for velocity
	Vector ro_tau_x(gr[1].x,gr[2].x,gr[3].x);				// ro*tau_x
	ro_tau_x = ro_tau_x + gr[1] - gr[0].x*U - U.x*gr[0];

	Vector ro_tau_y(gr[1].y,gr[2].y,gr[3].y);				// ro*tau_y
	ro_tau_y = ro_tau_y + gr[2] - gr[0].y*U - U.y*gr[0];
	
	Vector ro_tau_z(gr[1].z,gr[2].z,gr[3].z);				// ro*tau_z
	ro_tau_z = ro_tau_z + gr[3] - gr[0].z*U - U.z*gr[0];

	res[1] = (ny/a[0])*(((2/3)*ro_C*v_i) - ro_tau_x)*f.FaceNormal;
	res[2] = (ny/a[0])*(((2/3)*ro_C*v_j) - ro_tau_y)*f.FaceNormal;
	res[3] = (ny/a[0])*(((2/3)*ro_C*v_k) - ro_tau_z)*f.FaceNormal;
	//calc viscocity flux for energy
	Vector Grad_ro2_U2 = 2*(a[1]*gr[1] + a[2]*gr[2] + a[3]*gr[3]);
	Vector GradU2 = Grad_ro2_U2 - 2*a[0]*(U*U)*gr[0];
	GradU2 = (1/(a[0]*a[0]))*GradU2;
	Vector GradE = gr[4] - (a[4]/a[0])*gr[0];
	GradE = (1/a[0])*GradE;
	Vector p_e = 0.5*(model.gamma/Pr_com - 1)*GradU2 -(1/3)*(ro_C/a[0])*U - (model.gamma/(Pr_com))*GradE;
	res[4] = ny*p_e*f.FaceNormal;
	return res;
};

var_vector Method::VolumeSource(Cell &c) {
	var_vector Q(nv, 0);
	/*
	std::istringstream ss;	
	//Get all necessary parameters from list
	bool isExternalForce;
	std::string tmp = GetParam(ConfigParameters, "ExternalForce");
	ss.clear(); ss.str(tmp); ss >> isExternalForce;
	if (isExternalForce) {
		var_vector val = model.GetValue(c.GlobalIndex);
		std::string forceType = GetParam(ConfigParameters, "ExternalForceType");
		ss.clear(); ss.str(tmp); ss >> isExternalForce;
		if (forceType == "UniformGravity") {
			Vector g;
			Vector rov(val[1], val[2], val[3]);
			tmp = GetParam(ConfigParameters, "ExternalForceParams");
			ss.clear(); ss.str(tmp); ss >> g.x >> g.y >> g.z;
			Q[0] = 0;
			Q[1] = val[0] * g.x;
			Q[2] = val[0] * g.y;
			Q[3] = val[0] * g.z;
			Q[4] = rov*g;
		};
	};
	*/
	var_vector val = model.GetValue(c.GlobalIndex);
	Vector g(0,-1000,0);
	Vector rov(val[1], val[2], val[3]);

	Q[0] = 0;
	Q[1] = val[0] * g.x;
	Q[2] = val[0] * g.y;
	Q[3] = val[0] * g.z;
	Q[4] = rov*g;
	return Q;
};

double Method::Step(double time)
{
	//Define all initial steps	
	if (grad.size() != nc) grad.resize(nc);	
	//Calculate gradient in all cells if needed
	if ((SecondOrder) || (model.my != 0)) {
		#pragma omp parallel for
		for (int i = 0; i<nc; i++) {
			grad[i] = Grad(cells[i]);
		};
		if (P.NP > 1) {
			ExchangeGradValues();
		};
	};
	if (SecondOrder)
	{
		if (limit.size() != nc) limit.resize(nc);	
		#pragma omp parallel for
		for (int i = 0; i<nc; i++) {
			limit[i] = Limiter(cells[i]);
			//var_vector unity(nv,1);
			//limit[i] = unity;
		};
		if (P.NP > 1) {
			ExchangeLimiterValues();
		};
	};
	
	//Turbulent Block
	if ((time>=StartTurbulenceTime)&&(isBoundaryLayerActive)&&(model.my > 0)) {
		//#pragma omp parallel for
		for (std::map<int, CellBoundarySolution*>::iterator it = solutions.begin(); it!=solutions.end(); it++)
		{	
			CellBoundarySolution* sol = it->second;
			CalcP_grad(*sol, sol->up_face+1);
			//calc sigma and Flow_direction
			sol->Flow_direction = ChooseFlow_direction(*sol);

			if(CheckTurbulent(*sol))
            {
                sol->InitAproach();     //initial approaching
                sol->Solve(BoundaryC, BoundaryAlpha, BoundaryGamma, BoundaryRelativeTol, BoundaryMaxIter);//solver of 1D correlation equation 
				sol->U_inf = DefineUinf(*sol);	//calculate U_inf for connection between turbulent and laminar solution
				double sig = sol->U_inf*sol->ny_kinetic/(sol->h*sol->h);
				std::cout << "sigma_model=" << sol->sigma << " sigma_theory=" << sig << "\n";
				std::cin >> sig;
                sol->ComputeU_1();                //from knowed correlations compute U_averaged distribution
				sol->ComputeU_av();                 // and do it averaged for cells
				for(int i=0; i<sol->U_av.size(); i++)
				{
					int cell = CellInd(sol->upper_cell_ind[i]);
					//write the correct values for ro*V_turb
					values[cell][1] = sol->U_av[i]*values[cell][0]*sol->Flow_direction.x;
					values[cell][2] = sol->U_av[i]*values[cell][0]*sol->Flow_direction.y;
					values[cell][3] = sol->U_av[i]*values[cell][0]*sol->Flow_direction.z;
				};
            }else sol->TrivialSolution();            
		};
	};
	

	//Nullify fluxes			
	for (int i = 0; i<flux.size(); i++) flux[i] = var_vector(nv,0);		
	double time_step = 1e200;
	double max_sr = 0;	

	//Step through all faces to compute flux increment
	#pragma omp parallel for
	for (int i = 0; i<nf; i++) {
		double sr;
		//Solve reimann problem				
		std::pair<var_vector, double> res = Flux(faces[i]);
		flux[i] = res.first;
		sr = res.second;				

		if (model.my != 0) { 		
			
			flux[i] += ViscosityFlux(faces[i]);
			
			double h1 = model.GetCell(faces[i].FaceCell_1).CellVolume/faces[i].FaceSquare;
			h1 /= model.GetValue(faces[i].FaceCell_1)[0];
			double h2 = 0;
			if (!faces[i].isExternal) {
				h2 = model.GetCell(faces[i].FaceCell_2).CellVolume/faces[i].FaceSquare;
				h2 /= model.GetValue(faces[i].FaceCell_2)[0];
			};
			double h = h1>h2?h1:h2;
			sr += model.my*h; // TO DO viscosity constraint for time step
		};
		
		//Adjust time step
		if (P.IsLocal(faces[i].FaceCell_1)) {
			double tmp1 = cells[CellInd(faces[i].FaceCell_1)].CellVolume/sr;
			if (time_step>tmp1) time_step = tmp1;
		};
		if ((!faces[i].isExternal) && (P.IsLocal(faces[i].FaceCell_2))) {
			double tmp2 = cells[CellInd(faces[i].FaceCell_2)].CellVolume/sr;
			if (time_step>tmp2) time_step = tmp2;
		};		
	};	
	//Global minimum for time step definition
	time_step = model.CFL*P.GlobalMin(time_step);
	//Accumulate all fluxes to cells
	#pragma omp parallel for
	for (int i = 0; i<nf; i++) {		
		if (P.IsLocal(faces[i].FaceCell_1)) {
			var_vector flux_cell1 = time_step*(1.0/cells[CellInd(faces[i].FaceCell_1)].CellVolume)*flux[i];		
			values[CellInd(faces[i].FaceCell_1)] -= flux_cell1;
		};
		if ((!faces[i].isExternal) && (P.IsLocal(faces[i].FaceCell_2))) {
			var_vector flux_cell2 = time_step*(1.0/cells[CellInd(faces[i].FaceCell_2)].CellVolume)*flux[i];
			values[CellInd(faces[i].FaceCell_2)] += flux_cell2;
		};
	};
//	printf("Finished flux rank = %ld\n", P.rank);
	//Volume sources contribution	
	#pragma omp parallel for
	for (int i = 0 ; i<nc; i++) {
		Cell& c = model.GetCell(cells[i].GlobalIndex);		
		values[CellInd(c.GlobalIndex)] += time_step * CalcVolumeSources(c.GlobalIndex);
	};
	
	return time_step;
};

int Method::Run()
{	
	//Prepare storage
	flux.resize(nf, var_vector(nv, 0));
	residual.resize(values.size());
	//Open needed files
	std::ofstream resfile("residual.dat");
	double total_time = 0;
	if (P.NP > 1) {
		InitParallel();
	};
	if (P.IsMaster()) {
		printf("Started! NP = %ld\n", P.NP);
	};
	//calculation of faces array for all BoundaryCells for turbulent model
	if ((isBoundaryLayerActive) && (model.my > 0)) {
		//#pragma omp parallel for
		for (std::map<int, CellBoundarySolution*>::iterator it = solutions.begin(); it!=solutions.end(); it++)
		{		
			CellBoundarySolution* sol = it->second;
			for(int i=0; i<(sol->Max_faceNumber);i++)
			GetNextFace(*sol);
		};
		for (std::map<int, CellBoundarySolution*>::iterator it = solutions.begin(); it!=solutions.end(); it++)
		{		
			CellBoundarySolution* sol = it->second;
			double height=0;
			sol->h = 0.05;		//critical height
			sol->up_face=-1;
			for(int i=0; height<sol->h; i++)
			{
				height +=sol->cell_sizes[i];
				sol->up_face++;
			};
		};
	};

	double exec_time;
	for (int i = 0; i<Max_Iter; i++){
		if (P.NP > 1) {
			model.ExchangeValues();
		};		
		//printf("After exchange rank = %ld\n", P.rank);
		//Prepare residual
		if (i % OutputStep == 0) {
			for (int j = 0; j<values.size(); j++) residual[j] = values[j];
		};
		double dt = Step(total_time);
		//Calculate residual
		if (i % OutputStep == 0) {
			for (int j = 0; j<values.size(); j++) residual[j] = residual[j] - values[j];
		};
		total_time += dt;							
		if ((i % OutputStep == 0) && (!SuppressOutput) ) {
			//Calculate flow parameters
			double ro_av = 0;
			double ro_min = 1e200;
			double ro_max = -1e200;
			double Mach_max = 0;
			double volume = 0;
			for (int j=0; j<cells.size(); j++) {
				if (values[j][0] > ro_max) ro_max = values[j][0];
				if (values[j][0] < ro_min) ro_min = values[j][0];
				double M = Mach(values[j], model.gamma);
				if ( M > Mach_max) Mach_max = M;
				ro_av+= values[j][0] * cells[j].CellVolume;
				volume += cells[j].CellVolume;
			};
			ro_min = P.GlobalMin(ro_min);
			ro_max = P.GlobalMax(ro_max);
			ro_av = P.GlobalSum(ro_av);
			Mach_max = P.GlobalMax(Mach_max);
			volume = P.GlobalSum(volume);
			ro_av /= volume;
			//Temp
			double cold = 0;
			double hot = 0;
			for (int j=0; j<faces.size(); j++) {
				if ((faces[j].isExternal) && (faces[j].BCMarker == 2)) {					
				};
				if (faces[j].BCMarker == 3) hot += flux[j][4];
			};
			if (P.IsMaster()) {
				printf("Step: %ld , dt = %lg total_time = %lg\n",i,dt,total_time);	
				printf("ro_av : %lg\n", ro_av);
				printf("ro_min : %lg\n", ro_min);
				printf("ro_max : %lg\n", ro_max);
				printf("Mach_max : %lg\n", Mach_max);
			};
			//Calculate residual norm
			double norm2 = 0;
			for (int j = 0; j<residual.size(); j++) {
				norm2 += residual[j] * residual[j];
			};
			norm2 = P.GlobalSum(norm2);
			exec_time = P.GlobalTime();
			if (P.IsMaster()) printf("Execution time %lf minutes.\n", exec_time / 60);
			if (P.IsMaster()) printf("residual 2-norm = %lg\n", norm2);			
			resfile<<total_time<<" "<<norm2<<"\n";
			if (norm2 < ResidualTolerance) break;			
		};				
		

		double exec_time = P.GlobalTime();
		//if (P.IsMaster()) printf("Execution time %lf minutes.\n", exec_time / 60);

		if ((exec_time / 60.0) > MaxExecutionTime) {
			if (P.IsMaster()) printf("MaxExecutionTime reached.\n");
			return 0;		
		};
		if (total_time > t_max) return 0;
	};
	//Close files
	resfile.close();
	return 0;
};

//Load configuration
int Method::LoadConfig()
{	
	std::istringstream ss;	
	//Get all necessary parameters from list
	std::string tmp = GetParam(ConfigParameters, "MaxIter")[0];
	ss.clear(); ss.str(tmp); ss >> Max_Iter;
	tmp = GetParam(ConfigParameters, "MaxTime")[0];
	ss.clear(); ss.str(tmp); ss >> t_max;	
	tmp = GetParam(ConfigParameters, "MaxExecutionTime")[0];
	ss.clear(); ss.str(tmp); ss >> MaxExecutionTime;
	if (MaxExecutionTime == 0) MaxExecutionTime = 1e200;
	tmp = GetParam(ConfigParameters, "SuppressOutput")[0];
	ss.clear(); ss.str(tmp); ss >> SuppressOutput;
	tmp = GetParam(ConfigParameters, "OutputStep")[0];
	ss.clear(); ss.str(tmp); ss >> OutputStep;
	tmp = GetParam(ConfigParameters, "ResidualTolerance")[0];
	ss.clear(); ss.str(tmp); ss >> ResidualTolerance;
	tmp = GetParam(ConfigParameters, "SecondOrder")[0];
	ss.clear(); ss.str(tmp); ss >> SecondOrder;
	//Define volume sources
	InitVolumeSources();
	//Define boundary layer parameters
	InitBoundaryLayer();
	return 0;
};

//Read config and load parameters
void Method::InitVolumeSources() 
{
	std::istringstream ss;
	std::string type;
	VolumeSourcesTypes.clear();
	//Get all necessary parameters from list
	std::string tmp = GetParam(ConfigParameters, "VolumeSourceTypes")[0];		
	ss.clear(); ss.str(tmp); 
	isVolumeSourcesActive = false;
	while (!ss.eof()) {
		ss >> type;
		VolumeSourcesTypes[type] = true;
		isVolumeSourcesActive = true;
		type.clear();
	};
	//Init uniform gravity source
	if (VolumeSourcesTypes.find("UniformGravity") != VolumeSourcesTypes.end()) {
		tmp = GetParam(ConfigParameters, "UniformGravityVector")[0];
		ss.clear(); ss.str(tmp); ss >> g.x >> g.y >> g.z;
	} else {
		VolumeSourcesTypes["UniformGravity"] = false;
	};		
	return;
};

//Calculace volume sources term for given cell
var_vector Method::CalcVolumeSources(int GlobalIndex) {		
	var_vector Q(nv, 0);
	if (!isVolumeSourcesActive) return Q;
	if (VolumeSourcesTypes["UniformGravity"]) {
		var_vector val = model.GetValue(GlobalIndex);
		Vector rov(val[1], val[2], val[3]);
		Q[0] += 0;
		Q[1] += val[0] * g.x;
		Q[2] += val[0] * g.y;
		Q[3] += val[0] * g.z;
		Q[4] += rov*g;
	};
	return Q;
};


//Read config and init boundary layer info and parameters
void Method::InitBoundaryLayer() {
	int SolutionSize;
	std::istringstream ss;	
	//Get all necessary parameters from list
	std::string tmp = GetParam(ConfigParameters, "BoundaryLayer")[0];
	ss.clear(); ss.str(tmp); ss >> isBoundaryLayerActive;
	tmp = GetParam(ConfigParameters, "StartTurbulenceTime")[0];
	ss.clear(); ss.str(tmp); ss >> StartTurbulenceTime;
	if (isBoundaryLayerActive) {
		for (int i = 0; i<GetParam(ConfigParameters, "BoundaryLayerMarker").size(); i++) {
			int bc;
			tmp = GetParam(ConfigParameters, "BoundaryLayerMarker")[i];
			ss.clear(); ss.str(tmp); ss >> bc;
			BoundaryMarkers.push_back(bc);
		}
		tmp = GetParam(ConfigParameters, "BoundaryAlpha")[0];
		ss.clear(); ss.str(tmp); ss >> BoundaryAlpha;
		tmp = GetParam(ConfigParameters, "BoundaryC")[0];
		ss.clear(); ss.str(tmp); ss >> BoundaryC;
		tmp = GetParam(ConfigParameters, "BoundaryGamma")[0];
		ss.clear(); ss.str(tmp); ss >> BoundaryGamma;
		tmp = GetParam(ConfigParameters, "MaxBoundaryCells")[0];
		ss.clear(); ss.str(tmp); ss >> MaxBoundaryCells;
		tmp = GetParam(ConfigParameters, "BoundaryMaxIter")[0];
		ss.clear(); ss.str(tmp); ss >> BoundaryMaxIter;
		tmp = GetParam(ConfigParameters, "BoundaryRelativeTol")[0];
		ss.clear(); ss.str(tmp); ss >> BoundaryRelativeTol;		
		tmp = GetParam(ConfigParameters, "SolutionSize")[0];
		ss.clear(); ss.str(tmp); ss >> SolutionSize;
	}	

	//Init cell info
	for (int i = 0; i<faces.size(); i++) {
		bool isTurb = false;
		for (int j = 0; j<BoundaryMarkers.size(); j++) if ((faces[i].isExternal) && (faces[i].BCMarker == BoundaryMarkers[j])) {
			isTurb = true;
			break;
		};
		if (isTurb) {
			//New cell boundary solution object
			solutions[faces[i].GlobalIndex] = new CellBoundarySolution(SolutionSize, faces[i]);
			solutions[faces[i].GlobalIndex]->Max_faceNumber = MaxBoundaryCells;
		};
	};
	return;
};

void Method::SolveBoundaryLayer(Cell& c, Vector Uinf) {	
	return;
}

#endif